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Abstract 

We report first principles analysis of electron-phonon coupling in molecular devices under ex- 
ternal bias voltage and during current flow. Our theory and computational framework are based 
carrying out density functional theory within the Keldysh nonequilibrium Green's function formal- 
ism. We analyze which molecular vibrational modes are most relevant to charge transport under 
nonequilibrium conditions. For a molecular tunnel junction of a 1,4-benzenedithiolate molecule 
contacted by two leads, the low-lying modes of the vibration are found to be most important. As 
a function of bias voltage, the electron-phonon coupling strength can change drastically while the 
vibrational spectrum changes at a few percent level. 

PACS numbers: 81.07.Nb, 68.37.Ef, 72.10.-d, 73.63.-b 
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One of the most important questions concerning charge transport in molecular electronic 
devices is the role of electron-phonon (e-p) interaction. Here, "phonon" refers to quantized 
molecular vibrational modes which couple to various scattering states of the device. A 
typical molecular device has the Metal-Molecule-Metal (MMM) configuration schematically 
shown in Fig. Q where metal leads extend to far away and bias voltages can be applied so 
that a current flows through. The problem of predicting vibrational spectra and e-p coupling 
strength for such an open system during current flow, within self-consistent first principles 
including all atomic details of the molecule as well as the leads, is a serious theoretical 
challenge that has not been satisfactorily addressed. A particularly important problem is to 
understand which vibrational mode couples to which scattering state at what bias voltage 
. It is the purpose of this article to address this issue. 

Experimentally, single molecule vibrational spectra can be measured by inelastic tun- 
neling spectroscopy (IETS) H U, [J. Theoretically, various models have been applied to 
understand IETS and to investigate effects of e-p interaction based on tight binding atom- 
istic Hamiltonians [fj, Recently, Frederiksen et al. reported a first principles analysis 
of inelastic current due to e-p interactions in an Au chain, in which the relevant vibrations 
are along the chain length. In their theory [9|, the vibrational spectra was obtained using a 
plane-wave basis density functional theory (DFT) code in a cluster configuration at equilib- 
rium, and the dynamic matrix was evaluated using a finite differencing scheme. Transport 
properties were then obtained using the Transiesta package 10] with LCAO basis, and e-p 
scattering was included at the level of self-consistent Born approximation. 

In order to investigate voltage dependence of the e-p interaction in MMM devices, how- 
ever, quantized molecular vibrations and electrons need to be treated on equal footing at 



nonequilibrium. We accomplish this by carrying out D 
nonequilibrium Green's function (NEGF) formalism 



T atomic analysis within the Keldysh 



111 ]. In addition, we calculate the dy- 



namic matrix within the NEGF-DFT formalism ^ by evaluating analytical formula rather 
than numerical finite differencing ^: this is more general and more accurate so that all the 
phonon modes and e-p couplings can be obtained for complicated systems. For a molec- 
ular tunnel junction of a 1,4-benzenedithiolate (BDT) molecule contacted by two metallic 
electrodes (see Fig. [TJ, we found that the low- lying modes of the vibration are the most 
important for e-p coupling. As a function of bias voltage, the coupling strength can change 
drastically while the vibrational spectrum changes at a few percent level. 
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FIG. 1: Schematic plot of a Metal-BDT-Metal molecular tunneling junction (BDT = C6H4S2). 
The electrodes consist of repeated unit cells extending to z = ±00. The scattering region contains 
several layers of electrodes and the molecule. The vibrational box lies inside the scattering region. 

We start from the NEGF-DFT formalism documented in Ref. in which the density 
matrix p is calculated by NEGF G < , p ~ J dEG < (E). This way, we naturally take into 
account external bias voltage and open device transport boundary condition. After the 
Kohn-Sham Hamiltonian Hxs[p] of the device is obtained self-consistently, the total energy 
of the scattering region (see Fig. |TJ), as well as all transport properties of the device, can be 
obtained h|. Importantly, the NEGF-DFT formalism allows one to obtain Hrs and total 
energy E({Ri}, V&) as functions of external bias V&, (R, is the position of the i-th atom). 

The phonon (vibrational) eigenvectors e u and frequencies uo v [y is the mode index) are 
obtained by diagonalizing the dynamic matrix (Hessian matrix), 

U h y = V R ..V R . /J B({R i }, V h )/y/M^, (1) 

where M is the mass of an atom. Once (e u , uj u ) are obtained, the e-p interaction strength 



g u , denned by the e-p Hamiltonian 



12j, can be obtained from the standard expression: 



i 

Here, is a basis function for orbital p (p = s,p, d) of the j-th atom. The derivative of 

the KS Hamiltonian is carried out by fixing the atomic positions at their stationary locations 
R°. (R° is where the z-th atom feels no force ^^.) Although the Hessian matrix T~Lj } j> and 
the e-p coupling matrix g'ju-j' u 1 a PP ear to have the same form as those for equilibrium 
system, the derivatives of ^^^({Rj}, V&) and E({Ri},Vb) with respect to Rj propagate to 
derivatives of NEGF G < (E) that includes nonequilibrium physics [l3j . 



The e-p coupling is characterized by a dimensionless parameter 12j A e _ p which is con- 
tributed by all phonon modes, 
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Ae- P = J2 A "> X » = DOS(e F )^p^. (3) 
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Here, DOS(E) is the density of states of the scattering region and e F the Fermi energy of the 
leads. As we are interested in e-p coupling for quantum transport, the g v matrix is averaged 
over scattering states \l/ sc , which are obtained by the NEGF-DFT numerical package for 
any MMM device: 

(g»(E,E')) = {* sc (E)\g»\y sc (E')). (4) 

If both scattering states have the same energy E = E' = e F , we call such an e-p coupling 
the "elastic" one, A^. When E = e F and E' — e F ± fku^, we call it the "inelastic" A^. 
Finally, when the MMM device is under an external bias voltage we further average X u 
over the transport energy window (fi L , fi R ) where [i R / L are the electrochemical potentials of 
the right/left leads and \/i R — fi L \ — eVb- Hence, at nonequilibrium, we obtain 

J n L eVb huj v 

The inelastic coupling A^(Vb) is calculated by replacing g u (E,E) with g v (E, E ± fkv u ) in 
Eq. ©. 

The above theoretical formalism is implemented into our NEGF-DFT package McDCAL 
In numerical calculations, we further define a "vibrational box" inside the MMM device 
(see Fig. ^) which contains the atoms of interest. Typically, the vibrational box include the 
molecule and perhaps a few layers of the nearest lead atoms. The atomic indexes in the 
above formalism refer to those inside the vibrational box. 

In the following, we investigate general features of vibrational spectra and e-p coupling 
during nonequilibrium transport using the model MMM device of Fig. ^ i.e., a BDT 
molecular wire For each bias voltage, we iterate the KS Hamiltonian of the device 

to numerical convergence using the NEGF-DFT method [llj|; the atomic positions in the 
scattering region must also be relaxed for each applied bias. Afterward, the vibrational 
spectrum and the e-p coupling are obtained. As a check, we calculated u u of an isolated BDT 
using our NEGF-DFT formalism and obtained reasonable agreement, to within < 5 — 6 %, 
with experimental data collected by Raman spectroscopy and other methods fl5j . We also 



checked that the diagonal matrix elements of the e-p coupling are non-zero for modes having 
the A g symmetry and are zero for other modes, in agreement with selection rules from group 
theory 16]. 

When the BDT is placed between the leads (see Fig. [Q, new properties arise. First, 
several low-lying modes that do not exist for isolated BDT are found to play important roles. 
These modes include the center-of-mass and libration (CM(i) and LB(i), i = X,Y,Z) 
with energy hw u ~ 8 — 14 meV. Clearly, the presence of leads breaks the translational 
and rotational symmetries and produces these low-lying modes. Second, many vibrational 
frequencies are renormalized, to ~ 10 — 30 %, from that of the isolated molecule. This is 
especially true for modes with strong sulphur oscillations in the BDT. For large bias, Vb ~ 1 
V, we found that uj v changes up to 2 % while e u changes up to 5 % compared with the 
vibrational spectrum at Vb — 0. At Vb = 0, it turns out that all the modes can be classified 
by the same D 2 h point group as in the case of an isolated BDT; at 14 ^ 0, they can be 
classified by the C^v point group. 

Fig. Elplots the e-p coupling \ u versus ftw u at Vb = 0. For small bias voltages, Vb < 0.5 V, 
\ v does not change qualitatively. Most clearly shown is that some phonon modes give 
distinctly larger e-p coupling to scattering states than others (see Eq. (jlj)). Beside the 
expected in-plane A g (n) modes, modes of other symmetries are also responsible for the 
peaks in A„ (see the right insets of Fig. |2 for a few important modes). Notable are the 
in-plane modes Bi n (n), the center-of-mass mode CM(Z), the out-of-plane modes A n (n), 
-Big(l), and the libration LB(Z). Recall that for vibrational spectroscopy on free BDT 
such as Raman or infrared, there are always selection rules of modes |l5j. For a BDT 
device, however, our results suggest that no obvious selection rules are followed because 
many modes with very different symmetries manifest. Interestingly, among the low-lying 
modes with u v < 1000 cm -1 (huj u < 0.12 eV), the "breathing" modes A g (l) and A g {2) are 
not the most important ones for coupling to scattering states (A^ ~l-2x 10~ 4 ), although 
these modes are important for a free BDT. For the BDT device, the total e-p coupling (see 
Eq. ©) is Xf_ p « A~„ p « 7 x 10~ 3 at V b = 0. 

To understand why modes with symmetry other than A g can couple to scattering states, 
as shown by some of the peaks in Fig. |2l we consider Eq. (0J). For a free BDT, the 
electronic wave function in Eq. (@J) is a molecular orbital, and for each orbital one obtains 
a value (g u )- Hence for a free BDT the relevant e-p coupling is a diagonal matrix in orbital 
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FIG. 2: Dimensionless e-p coupling constants Xfj (circles), A~ (triangles), A+ (diamonds) as a 
function of vibrational frequency uj u at V& = 0. The lines are guide to the eye. The modes are 
classified using the D^h point group. The left inset shows projection P sc ,rmo of a scattering 
state at Fermi energy (shifted to zero) onto molecular orbitals (RMO) of the BDT. Letters H and 
L mean HOMO and LUMO. In the right side insets, we show eigenvectors of several important 
eigenmodes of the BDT at small biases, namely, the libration LB(Z) (top), ^4 g (4) (middle), and 
jBiu(4) (bottom). 

space: only those vibrational modes with A, symmetry give nonzero va.nes to these diagonal 
matrix elements |16(. For transport, however, the wave function appearing in Eq. (@J) is a 
scattering state which is roughly a linear combination of many molecular orbitals. Therefore 
it is possible to have off-diagonal matrix elements in the coupling matrix so that modes with 
symmetries other than A g can also contribute. This can be substantiated as follows. We 
project scattering states $? SC (E) onto "renormalized molecular orbitals" (RMO) of the BDT 
in the MMM device The RMO's are obtained by diagonalizing the Hamiltonian sub- 
matrix that corresponds to the BDT molecule, and this sub-matrix is a part of the total 
KS Hamiltonian of the entire MMM device hj. Note that RMO's can be different from 
the original molecular orbitals of an isolated BDT due to charge transfer from the leads to 
the molecule and external bias potentials. The projection is characterized by the quantity 
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FIG. 3: Dimensionless e-p coupling constants (circles) and A~ (triangles) at 14 = 1 V. Although 
in principle one should use the C 2 v point group to label the modes, for comparison with Fig. |2] we 
keep the D 2 h labels here. The left insets show projection P SCi rmo for two scattering states taken 
at E — e F = 0.4 eV. In the right side insets, we show eigenvectors of several important modes at 
V b ~ 1 V, namely, CM{Y) (top), CM{Z) (middle), and A g (l) (bottom). 

Psc,rmo = |(^sc(£f) I RMO)\ 2 plotted in the left inset of Fig. El as a function of energy. 
We found that all scattering states \l/ sc near the Fermi energy are contributed by the same 
several dominant RMO's at low bias. Therefore one can well consider that \l/ sc is a linear 
combination of these few RMO's and g u of Eq. (J3J) is contributed mostly by them, i.e., ( g v ) 
is contributed by a quantity g v aa , = (RMO a \g v '\RMO a >). Hence, for transport problems, 
the off-diagonal contributions (when a ^ a') can be as important as the diagonal ones 
a — a'). Furthermore, when bias V\, is increased, molecular orbitals in the MMM device 
1t| become less symmetric so that vibrational modes different from the A g symmetry can 
even contribute to the diagonal matrix elements of the e-p coupling. For example, at 14 = 1 
V, the contribution of orbitals HOMO — n with n — 1, 3 to the scattering states is found to 
give rise to non-zero e-p coupling for vibrational modes with both the A g and Bi u symmetries. 
These results allow us to conclude that e-p coupling in MMM devices during current flow 
(Vb ^ 0) is much more complicated than that for free molecules. 
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A most important finding is that at large bias voltage, V& ~ 1 V, the e-p coupling 
strength changes drastically although the vibrational spectrum is only changed by a few 
percent as mentioned above. Fig. |3] plots the coupling strength at V& = 1 V. We observe 
that contributions to e-p coupling are now dominated by a few low-lying modes such as the 
center-of-mass modes; the total coupling is a factor of five greater than that at low bias (see 
Fig. |2J), and coupling due to individual low- lying modes is also much larger. By projecting 
different scattering states ^> SC (E) with the same energy E onto RMO's as presented above, 
we found that two patterns of P sc ,rmo occur, as shown in the left insets of Fig. El This 
indicates that scattering states inside the transport energy window are different from those 
at low bias. In particular, the pattern of P sc ,rmo m the upper panel is similar to that in 
Fig. |21 but the new (lower) pattern of P SCt rmo comes from lower HOMO — n and higher 
LUMO + n RMO's. This leads to a different behavior of e-p couplings at large biases for the 
BDT device. For example, our calculations reveal that the peak labeled CM(Y) in \ v {u v ) 
(see Fig. comes from particular off-diagonal matrix elements, (L + 2 | g v \ H — 1) and 
(L + 4 | g u | H — 1). The peak labeled CM(Z), on the other hand, is found to come from 
diagonal matrix elements, e.g., (H — 1 1 g v \ H — 1). These findings also correlate well with 
the peaks of P sc ,rmo i n the left insets of Fig. El The total A c _ p at VJ, — 1 V is found to 
be ~ 0.04. This enhancement by roughly a factor of five from that of VJ, ~ is due to the 
center-of-mass modes shown in Fig. El that can be confirmed by computing A e _ p without 
counting these modes. 

Why bias voltage can change e-p coupling so drastically? We found that the reason is 
mainly due to contribution of different scattering states. In the inset of Fig. EJ we plot 
transmission coefficient T(E, Vj>) vs. electron energy E for several values of V&. Most clearly 
shown is that V& shifts the transmission features toward the transport window. In particular, 
a sharp peak (at E ~ —0.3 eV) is shifted up-wards in energy with the increase of V&. When 
Vb > 0.5 V, the "tail" of this sharp peak starts entering the transport window (between 
\i L = and fj JR = eVb, see also Eq. (JHJ)). If V& > 0.75 V, this peak enters into the transport 
window completely. When this happens, the e-p coupling changes drastically and the new 
pattern appears in the projection P sc ,rmo as discussed above. Fig. |U plots the e-p coupling 
strength A„ for several vibrational modes v and the total A e _ p vs. V&. The curves give a 
clear indication that the e-p coupling is roughly a constant at small biases, but can change 
nonlinearly as the applied bias voltage is varied. Such a change can have deep implications 
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FIG. 4: Dimensionless e-p coupling constants \f_ p (circles) and , for v = CM{Y) (squares), 
CM(Z) (diamonds), and A g (l) (triangles), as a function of V b . Inset: transmission function 
T(E,V b ) at V b = 0.5 V (full line), V b = 0.75 V (dashed line), and V b = 1 V (dash-dot line). The 
transport energy window (n L , /x R ) lies in the positive part of the energy axis, from E = to 
E=\eV b \. 



to local heating in the device during nonequilibrium charge transport 




In summary, the entire relevant vibrational spectrum of a molecule device can be obtained 
at non-zero bias within the NEGF-DFT formalism where both vibrational and electronic 
properties are calculated at equal footing. For a 1,4-BDT molecular device studied here, low- 
lying vibrational modes play an important role in contributing to the e-p coupling strength. 
The coupling strength changes drastically as bias voltage is increased due to participations 
of new scattering states. For the BDT device, at large bias of V b ~ 1 V, it is the center- 
of-mass modes which denominate the e-p coupling, while for small bias V b < 0.5 V, many 
modes of different symmetries contribute. The vibrational spectrum also depends on bias, 
but for the BDT device this dependence is at a few percent level. 
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